R Markdown
setwd("/Users/macbook/Desktop/Bio-Research\ /PS050_cellranger_count_outs/raw_feature_bc_matrix")
EpCAM.data <- Read10X(data.dir = "/Users/macbook/Desktop/Bio-Research\ /PS050_cellranger_count_outs/filtered_feature_bc_matrix/")
# Create Seurat object
EpCAM <- CreateSeuratObject(counts = EpCAM.data, project = "mSG_EpCAM", save.SNN = TRUE)
EpCAM
## An object of class Seurat
## 32285 features across 8277 samples within 1 assay
## Active assay: RNA (32285 features, 0 variable features)
# Add number of genes per UMI for each cell to metadata
# EpCAM$log10GenesPerUMI <- log10(EpCAM$nFeature_RNA) - log10(EpCAM$nCount_RNA) ?
EpCAM$log10GenesPerUMI <- log10(EpCAM$nFeature_RNA) / log10(EpCAM$nCount_RNA)
# Compute percent mito ratio
# What is the meaning of `pattern = "^mt-" ` here?
EpCAM$mitoRatio <- PercentageFeatureSet(object = EpCAM, pattern = "^mt-")
EpCAM$mitoRatio <- EpCAM@meta.data$mitoRatio / 100
# Create metadata dataframe
metadata <- EpCAM@meta.data
# Add cell IDs to metadata
metadata$cells <- rownames(metadata)
# Rename columns
metadata <- metadata %>%
dplyr::rename(seq_folder = orig.ident,
nUMI = nCount_RNA,
nGene = nFeature_RNA)
# Create sample column
metadata$sample <- "mSG"
# Create metadata dataframe
metadata <- EpCAM@meta.data
# Add cell IDs to metadata
metadata$cells <- rownames(metadata)
# Rename columns
metadata <- metadata %>%
dplyr::rename(seq_folder = orig.ident,
nUMI = nCount_RNA,
nGene = nFeature_RNA)
# Create sample column
metadata$sample <- "mSG"
# Add metadata back to Seurat object
EpCAM@meta.data <- metadata
# Cell-level FILTERING
set.seed (12)
filtered_seurat <- subset(x = EpCAM,
subset= (nUMI >= 1) &
(nUMI <=20000) &
(nGene >= 100) &
(log10GenesPerUMI > 0.5) &
(log10GenesPerUMI < 0.9) &
(mitoRatio < 0.20))
# Output a logical vector for every gene on whether the more than zero counts per cell
# Extract counts
counts <- GetAssayData(object = filtered_seurat, slot = "counts")
# Save filtered subset to new metadata
metadata_clean <- filtered_seurat@meta.data
filtered_seurat[["percent.mt"]] <- PercentageFeatureSet(filtered_seurat, pattern = "^mt-")
plot1 <- FeatureScatter(filtered_seurat, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(filtered_seurat, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2+ plot_layout(ncol = 1)
# Normalize the counts
set.seed(12)
filtered_seurat_norm <- SCTransform(filtered_seurat, vars.to.regress = "mitoRatio")
## Calculating cell attributes from input UMI matrix: log_umi
## Variance stabilizing transformation of count matrix of size 14745 by 2806
## Model formula is y ~ log_umi
## Get Negative Binomial regression parameters per gene
## Using 2000 genes, 2806 cells
##
|
| | 0%
|
|================== | 25%
|
|=================================== | 50%
|
|==================================================== | 75%
|
|======================================================================| 100%
## Found 108 outliers - those will be ignored in fitting/regularization step
## Second step: Get residuals using fitted parameters for 14745 genes
##
|
| | 0%
|
|== | 3%
|
|===== | 7%
|
|======= | 10%
|
|========= | 13%
|
|============ | 17%
|
|============== | 20%
|
|================ | 23%
|
|=================== | 27%
|
|===================== | 30%
|
|======================= | 33%
|
|========================== | 37%
|
|============================ | 40%
|
|============================== | 43%
|
|================================= | 47%
|
|=================================== | 50%
|
|===================================== | 53%
|
|======================================== | 57%
|
|========================================== | 60%
|
|============================================ | 63%
|
|=============================================== | 67%
|
|================================================= | 70%
|
|=================================================== | 73%
|
|====================================================== | 77%
|
|======================================================== | 80%
|
|========================================================== | 83%
|
|============================================================= | 87%
|
|=============================================================== | 90%
|
|================================================================= | 93%
|
|==================================================================== | 97%
|
|======================================================================| 100%
## Computing corrected count matrix for 14745 genes
##
|
| | 0%
|
|== | 3%
|
|===== | 7%
|
|======= | 10%
|
|========= | 13%
|
|============ | 17%
|
|============== | 20%
|
|================ | 23%
|
|=================== | 27%
|
|===================== | 30%
|
|======================= | 33%
|
|========================== | 37%
|
|============================ | 40%
|
|============================== | 43%
|
|================================= | 47%
|
|=================================== | 50%
|
|===================================== | 53%
|
|======================================== | 57%
|
|========================================== | 60%
|
|============================================ | 63%
|
|=============================================== | 67%
|
|================================================= | 70%
|
|=================================================== | 73%
|
|====================================================== | 77%
|
|======================================================== | 80%
|
|========================================================== | 83%
|
|============================================================= | 87%
|
|=============================================================== | 90%
|
|================================================================= | 93%
|
|==================================================================== | 97%
|
|======================================================================| 100%
## Calculating gene attributes
## Wall clock passed: Time difference of 37.43277 secs
## Determine variable features
## Place corrected count matrix in counts slot
## Regressing out mitoRatio
## Centering data matrix
## Set default assay to SCT
# Check which assays are stored in objects
filtered_seurat_norm@assays
## $RNA
## Assay data with 32285 features for 2806 cells
## First 10 features:
## Xkr4, Gm1992, Gm19938, Gm37381, Rp1, Sox17, Gm37587, Gm37323, Mrpl15,
## Lypla1
##
## $SCT
## SCTAssay data with 14745 features for 2806 cells, and 1 SCTModel(s)
## Top 10 variable features:
## Klk1, Wfdc18, Smgc, Klk1b26, Pip, Scgb1b27, Scgb2b27, Car6, Klk1b11,
## Mucl2
# Identification of highly variable features (feature selection)
filtered_seurat_norm <- FindVariableFeatures(filtered_seurat_norm, selection.method = "vst", nfeatures = 2000)
# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(filtered_seurat_norm), 10)
par(plt = c(0.1, 0.9, 0.1, 0.9))
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(filtered_seurat_norm)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
## When using repel, set xnudge and ynudge to 0 for optimal results
plot1 + plot2+ plot_layout(ncol = 1)
# Run PCA
set.seed (12)
filtered_seurat_norm_PCA <- RunPCA(object = filtered_seurat_norm)
## PC_ 1
## Positive: S100a6, Fxyd2, Nupr1, Actg1, Mt1, Anxa1, Crisp3, Cldn4, Stap1, Ubc
## Ifi27l2a, Actb, Krt18, Tnfrsf12a, Ifrd1, Atf3, Krt8, Ifitm3, Mt2, Ascl3
## Maff, Ubb, S100a10, Ier3, Gadd45a, Ier5, Cebpd, Ly6d, Krt7, Ldhb
## Negative: Pip, Wfdc12, Smgc, A630073D07Rik, Scgb2b27, Scgb1b27, Aqp5, Mucl2, Prol1, Car6
## Scgb2b26, Cldn10, Prlr, Phyh, Lars2, Lman1l, Gm44410, St3gal4, Serpinb11, Tll1
## Gm26917, Agt, Col11a1, Lrrc26, Fxyd3, Gdpd1, Hp, Rexo2, H2afj, Fam25c
## PC_ 2
## Positive: Klk1, Klk1b11, Klk1b26, Klk1b5, Klk1b4, Klk1b9, Bglap3, Klk1b3, Klk1b16, Ly6a
## Lgals3, Fosb, Scgb1c1, Crisp3, Klk1b22, Cyp2f2, Klk1b21, Timp3, Ces1d, Pcp4l1
## Klk1b27, Krt8, Foxq1, Arg1, Serpinb1a, Cldn4, Ngf, Egf, Fos, G0s2
## Negative: Stap1, S100a6, Ifi27l2a, Slc12a2, Ascl3, Hepacam2, Nupr1, Slc4a11, Ifitm3, Aldh1a3
## Pam, Foxi1, Serpinb9, Itih5, Pip, Rcan2, Ncald, Asgr1, Tspan1, Smgc
## Ldhb, Cebpd, Cpxm2, Vps36, Itpr2, Nrip3, Pdlim3, Aldh1a1, Rbpms, A630073D07Rik
## PC_ 3
## Positive: Wfdc18, Esp4, Car2, Tmem176b, Tmem176a, Ccnd1, Gjb2, Cyp2f2, Ctsh, Mgst1
## Hp, Aox3, BC006965, Ano1, Esp18, Elf5, Pir, Barx2, Slc13a2, Ehf
## Fth1, Aldoc, Palmd, Ces1d, Cdo1, Cyba, Cebpb, Taldo1, Lmo4, Nkd2
## Negative: Klk1, Crisp3, Klk1b26, Fxyd2, Klk1b11, Klk1b5, Klk1b9, Klk1b4, Smgc, Pip
## Mt1, Klk1b3, A630073D07Rik, Phyh, Klk1b22, Klk1b16, Aqp5, Mt2, Scgb1c1, Wfdc12
## Klk1b21, Cldn10, Klk1b27, Stap1, S100a6, Ngf, Wfdc2, Serpinb1a, Lman1l, Gm44410
## PC_ 4
## Positive: Crisp3, Fxyd2, Wfdc18, Ifi27l2a, Bglap3, Nupr1, Car2, Krt18, Anxa1, Stap1
## Esp4, Cyp2f2, Krt8, Ascl3, Slc12a2, Klk1b11, Klk1b26, Ldhb, Cldn4, Hepacam2
## Cryab, Gjb2, Ckmt1, Serpinb6b, Cited4, Cxcl17, Ncald, Actg1, Tmsb4x, Serpinb9
## Negative: Dcn, Sparc, Smoc2, Apoe, Rbp1, Igfbp4, Mgp, Col4a1, Vim, Serpinh1
## Lgals1, Htra1, Igfbp7, Cd34, Ebf1, Mdk, Rnase4, Hspg2, Cxcl12, Emp3
## Ppic, Aebp1, Plpp3, Efemp1, Lrp1, AW112010, C1s1, Lamb1, Fbln2, Meg3
## PC_ 5
## Positive: Actg1, Ifrd1, Krt8, Arid5b, Krt18, Actb, Tnfrsf12a, Cldn4, Pip, S100a10
## Ly6d, Sfn, Smgc, Tuba1c, Maff, Fosl1, Lmna, Fhl2, H3f3b, Tacstd2
## Ndrg1, Mast4, Cdkn1a, A630073D07Rik, Cldn10, Sprr1a, Ptn, Hspb8, Wfdc12, Phyh
## Negative: Crisp3, Ifi27l2a, Klk1b5, Fxyd2, Klk1b4, G0s2, Klk1b26, Klk1, Klk1b9, Klk1b22
## Fos, Klk1b16, Klk1b3, Stap1, Ascl3, Scgb1c1, Esp4, Klk1b21, Bglap3, Klk1b11
## Itih5, Tmem176b, Tmem176a, Klk1b27, Ngf, Cd81, Car2, Ldhb, Hepacam2, Hspa1a
# Plot PCA
PCAPlot(filtered_seurat_norm_PCA)
# Run UMAP
set.seed(12)
filtered_seurat_norm_UMAP <- RunUMAP(filtered_seurat_norm_PCA,
dims = 1:15,
reduction = "pca")
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 16:23:16 UMAP embedding parameters a = 0.9922 b = 1.112
## 16:23:16 Read 2806 rows and found 15 numeric columns
## 16:23:16 Using Annoy for neighbor search, n_neighbors = 30
## 16:23:16 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 16:23:17 Writing NN index file to temp file /var/folders/v2/mpyxd00n14zdpng47ksgc4z80000gn/T//RtmpNUNMgi/file1408a1f654025
## 16:23:17 Searching Annoy index using 1 thread, search_k = 3000
## 16:23:17 Annoy recall = 100%
## 16:23:18 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 16:23:19 Initializing from normalized Laplacian + noise (using irlba)
## 16:23:19 Commencing optimization for 500 epochs, with 112676 positive edges
## 16:23:23 Optimization finished
# Plot UMAP
DimPlot(filtered_seurat_norm_UMAP)
# Run t-SNE
filtered_seurat_norm_tSNE <- RunTSNE(
filtered_seurat_norm_PCA, reduction = "pca",
assay = NULL,
seed.use = 12,
tsne.method = "Rtsne",
dim.embed = 2,
dims = 1:20,
reduction.key = "tSNE_")
# Plot t-SNE
TSNEPlot(filtered_seurat_norm_tSNE)
# CLUSTERING CELLS BASED ON TOP PCs (METAGENES)
# Identify significant PCs
DimHeatmap(filtered_seurat_norm_tSNE,
dims = 1:12,
cells = 500,
balanced = TRUE)
# Printing out the most variable genes driving PCs
print(x = filtered_seurat_norm_tSNE[["pca"]],
dims = 1:12,
nfeatures = 15)
## PC_ 1
## Positive: S100a6, Fxyd2, Nupr1, Actg1, Mt1, Anxa1, Crisp3, Cldn4, Stap1, Ubc
## Ifi27l2a, Actb, Krt18, Tnfrsf12a, Ifrd1
## Negative: Pip, Wfdc12, Smgc, A630073D07Rik, Scgb2b27, Scgb1b27, Aqp5, Mucl2, Prol1, Car6
## Scgb2b26, Cldn10, Prlr, Phyh, Lars2
## PC_ 2
## Positive: Klk1, Klk1b11, Klk1b26, Klk1b5, Klk1b4, Klk1b9, Bglap3, Klk1b3, Klk1b16, Ly6a
## Lgals3, Fosb, Scgb1c1, Crisp3, Klk1b22
## Negative: Stap1, S100a6, Ifi27l2a, Slc12a2, Ascl3, Hepacam2, Nupr1, Slc4a11, Ifitm3, Aldh1a3
## Pam, Foxi1, Serpinb9, Itih5, Pip
## PC_ 3
## Positive: Wfdc18, Esp4, Car2, Tmem176b, Tmem176a, Ccnd1, Gjb2, Cyp2f2, Ctsh, Mgst1
## Hp, Aox3, BC006965, Ano1, Esp18
## Negative: Klk1, Crisp3, Klk1b26, Fxyd2, Klk1b11, Klk1b5, Klk1b9, Klk1b4, Smgc, Pip
## Mt1, Klk1b3, A630073D07Rik, Phyh, Klk1b22
## PC_ 4
## Positive: Crisp3, Fxyd2, Wfdc18, Ifi27l2a, Bglap3, Nupr1, Car2, Krt18, Anxa1, Stap1
## Esp4, Cyp2f2, Krt8, Ascl3, Slc12a2
## Negative: Dcn, Sparc, Smoc2, Apoe, Rbp1, Igfbp4, Mgp, Col4a1, Vim, Serpinh1
## Lgals1, Htra1, Igfbp7, Cd34, Ebf1
## PC_ 5
## Positive: Actg1, Ifrd1, Krt8, Arid5b, Krt18, Actb, Tnfrsf12a, Cldn4, Pip, S100a10
## Ly6d, Sfn, Smgc, Tuba1c, Maff
## Negative: Crisp3, Ifi27l2a, Klk1b5, Fxyd2, Klk1b4, G0s2, Klk1b26, Klk1, Klk1b9, Klk1b22
## Fos, Klk1b16, Klk1b3, Stap1, Ascl3
## PC_ 6
## Positive: Fos, Zfp36, Atf3, Junb, Smgc, Jun, Btg2, Egr1, Ier2, Irf1
## Gadd45g, H3f3b, Fosb, Ier3, Gadd45b
## Negative: Mucl2, Prol1, Scgb1b27, Car6, Scgb2b26, Scgb2b27, S100a6, Ly6d, S100a10, Ptn
## Agt, Dcn, Pigr, Ndrg1, F3
## PC_ 7
## Positive: Smgc, Phyh, Rps8, Fxyd3, Cyp2f2, Rpl13, A630073D07Rik, Cldn10, Prlr, Hp
## Eef1a1, Gm44410, Rplp0, Rpl32, Serpinb11
## Negative: Mucl2, Scgb1b27, Prol1, Car6, Scgb2b27, Scgb2b26, Atf3, Fos, Agt, Egr1
## Fosb, Gm26917, Jun, Pigr, Lpo
## PC_ 8
## Positive: Ifrd1, Krt8, Krt18, Anxa1, Cldn4, Bglap3, Actg1, Clu, Krt23, Pip
## Arid5b, Lgals3, Actb, Wfdc12, Ubc
## Negative: Fos, Ly6d, Jun, Ptn, Krt14, Dusp1, Hspa1a, Ucma, Ccn1, Krt5
## Lgals7, Krt17, Dcn, Ier2, Smoc2
## PC_ 9
## Positive: Isg15, Iigp1, Ifit1, Ifitm3, Ly6a, Pecam1, Srgn, Fabp4, Ecscr, Tm4sf1
## Rsad2, Emcn, Rasip1, Cd93, Flt1
## Negative: Dcn, Fos, S100a6, Egr1, Atf3, Nupr1, Krt18, Smoc2, Fosb, Mt1
## Jun, Rbp1, Gsn, Esp18, Klf6
## PC_ 10
## Positive: Klk1b22, Klk1b4, Klk1b5, Klk1b21, Egf, Ngf, Klk1b26, Klk1b3, Gm42418, Gm26917
## Wfdc18, Scgb1c1, Cracr2a, Esp4, Ier5
## Negative: Fxyd2, Cat, Ly6a, Bglap3, Fos, Arg1, Ces1d, Gstm1, Sord, Duoxa2
## Tacstd2, Krt7, Cyp2f2, Duox2, Wfdc2
## PC_ 11
## Positive: Gm26917, Gm42418, Malat1, Lars2, Cracr2a, Ppp1r10, Hspa1b, Neat1, Mafb, Hexim1
## Klk1, Jun, Clec2d, 2900060B14Rik, Hspa1a
## Negative: Atf3, Klk1b4, Pip, Mt1, Agt, Wfdc12, Bglap3, Klk1b5, Irf1, Ifitm3
## Zfp36, Crisp3, Nfkbia, Ier3, Cxcl10
## PC_ 12
## Positive: Hspa1a, Hspb1, Actg1, Hspa1b, Esp18, Actb, Ly6d, Hspe1, Timp3, Ubc
## Nupr1, Krt8, Hspa8, Fxyd2, Hsp90aa1
## Negative: Bglap3, Cxcl10, Gm26917, Irf1, Isg15, Ifi27l2a, Ifitm3, Ifit1, Irf7, Nfkbia
## Gm42418, Hsd11b1, Oasl2, Lars2, Duox2
VizDimLoadings(filtered_seurat_norm_tSNE, dims = 1, reduction = "pca")
VizDimLoadings(filtered_seurat_norm_tSNE, dims = 2, reduction = "pca")
VizDimLoadings(filtered_seurat_norm_tSNE, dims = 3, reduction = "pca")
# Elbow plot: threshold for identifying the majority of the variation. However, this method can be quite subjective.
ElbowPlot(object = filtered_seurat_norm_UMAP,
ndims = 40)
# CLUSTER THE CELLS
# Determine the K-nearest neighbor graph
set.seed (12)
filtered_seurat_norm_UMAP <- FindNeighbors(object = filtered_seurat_norm_UMAP, reduction = "pca", dims = 1:15, verbose = TRUE, graph.name = "mSG")
## Computing nearest neighbor graph
## Computing SNN
## Only one graph name supplied, storing nearest-neighbor graph only
# Determine the clusters with Leiden algorithm
set.seed (12)
filtered_seurat_norm_UMAP_0.5 <- FindClusters(object = filtered_seurat_norm_UMAP, verbose = TRUE, algorithm = 4, resolution = 0.5, graph.name = "mSG")
set.seed(12)
distance_matrix <- dist(Embeddings(filtered_seurat_norm_UMAP_0.5[['umap']])[, 1:2])
#Isolate cluster name
clusters <- filtered_seurat_norm_UMAP_0.5@active.ident
#Compute silhouette score for each cluster
silhouette <- silhouette(as.numeric(clusters), dist = distance_matrix)
filtered_seurat_norm_UMAP_0.5@meta.data$silhouette_score <- silhouette[,3]
#Plot
fviz_silhouette(silhouette, label = FALSE, print.summary = TRUE)
## cluster size ave.sil.width
## 1 1 672 0.75
## 2 2 474 0.18
## 3 3 370 0.66
## 4 4 364 0.81
## 5 5 255 0.40
## 6 6 221 0.53
## 7 7 168 0.61
## 8 8 151 -0.16
## 9 9 131 0.63
# Determine metrics to plot present in seurat_integrated@meta.data
metrics <- c("nUMI", "nGene", "mitoRatio")
FeaturePlot(filtered_seurat_norm_UMAP_0.5,
reduction = "umap",
features = metrics,
pt.size = 0.4,
order = TRUE,
min.cutoff = 'q10',
label = TRUE)
# Basal - Cluster 8
FeaturePlot(filtered_seurat_norm_UMAP_0.5_RNA,
reduction = "umap",
features = c("Krt15", "Krt14"),
order = TRUE,
min.cutoff = 'q10',
label = TRUE,
repel = TRUE)
# ACINAR 1 - Cluster 4
FeaturePlot(filtered_seurat_norm_UMAP_0.5_RNA,
reduction = "umap",
features = c("Pip", # for 1 and 2
"Scgb2b27", "Lpo", "Car6"), # for 1
order = TRUE,
min.cutoff = 'q10',
label = TRUE,
repel = TRUE)
# ACINAR 2
FeaturePlot(filtered_seurat_norm_UMAP_0.5_RNA,
reduction = "umap",
features = c("Pip", # for 1 and 2
"Prlr", "Smgc"), # for 2
order = TRUE,
min.cutoff = 'q10',
label = TRUE,
repel = TRUE)
# Ductal 1
FeaturePlot(filtered_seurat_norm_UMAP_0.5_RNA,
reduction = "umap",
features = c("Klk1",
"Ngf"),
order = TRUE,
min.cutoff = 'q10',
label = TRUE,
repel = TRUE)
# Ductal 2
FeaturePlot(filtered_seurat_norm_UMAP_0.5_RNA,
reduction = "umap",
features = c("Foxi1",
"Ascl3"),
order = TRUE,
min.cutoff = 'q10',
label = TRUE,
repel = TRUE)
# Ductal 3
FeaturePlot(filtered_seurat_norm_UMAP_0.5_RNA,
reduction = "umap",
features = c("Slc9a3", "Clic6"),
order = TRUE,
min.cutoff = 'q10',
label = TRUE,
repel = TRUE)
